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ABSTRACT 

We use time-dependent, one-dimensional disc models to investigate the evolution of 
protostellar discs that form through the collapse of molecular cloud cores and in which 
the primary transport mechanism is self-gravity. We assume that these discs settle 
into a state of thermal equilibrium with Q = 2 and that the strength of the angular 
momentum transport is set by the cooling rate of the disc. The results suggest that 
these discs will attain a quasi-steady state that persists for a number of free-fall times 
and in which most of the mass within 100 au is located inside 10 — 20 au. This pile-up of 
mass in the inner disc could result in temperatures that are high enough for the growth 
of MHD turbulence which could rapidly drain the inner disc and lead to FU Orionis- 
like outbursts. In all our simulations, the inner regions of the discs (r < 40 au) were 
stable against fragmentation, while fragmentation was possible in the outer regions 
(r > 40 au) of discs that formed from cores that had enough initial angular momentum 
to deposit sufficient mass in these outer regions. The large amounts of mass in these 
outer regions, however, suggests that fragmentation will lead to the formation of sub- 
stellar and stellar mass companions, rather than planetary mass objects. Although 
mass accretion rates were largely consistent with observations, the large disc masses 
suggest that an additional transport mechanism (such as MRI occuring in the upper 
layers of the disc) must operate in order to drain the remaining disc material within 
observed disc lifetimes. 

Key words: stars: formation — stars: pre-main-sequence — circumstellar matter — 
planetary systems: protoplanetary discs — planetary systems: formation 



1 INTRODUCTION 

The formation of low-mass stars occ urs through the collapse 
of col d, dense molecular cloud cores (iTerebev. Shu fc Cassenl 
1 1984I ) . These cores, however, generally contain amounts of 
angular momentum far in excess of the rotat ional angular 
momentum of a single star l|Caselli et al ]|2002). Most of the 
mass must therefore first pass through a circumstellar disc 
before accreting onto the central star. One of the open prob- 
lems in star and planet formation is what mechanism acts 
to transport the angular momentum outwards, allowing this 
accretion to take place. 

It is now generally accepted that in most astro- 
physical discs, angular momentum transport is driven 
by magnetohydrodynamic (MHD) turbu lence initiated by 
the magneto-rotational instability (MRI) (|f3albus fc Hawlevl 
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ll99ll : |Papal oizou fc NelsorJliool ). During the earliest stages 
of star formation protoplanetary discs are, however, cold 
and dense and probably do not have even the very small 
degree of ioniza t ion n eeded to sustain MHD turbulence 
l|Blaes fc Balbusl I1994T ). At these early times, the disc 
is, however, likely to be massive and disc self-gravity 
may then provide an alternate and possibly dominant 
transport mechanism through the growth of the grav- 
itational instability (iToomrei 1 1964 iLin fc Pringie] 1 19871 : 
lLaughlin fc Bodenheimerll 19941 ) . 

Although it has been suggested that gravitationally 
unstable discs may fragme nt to form gas g iant planets 
or substellar companions (Boss 19981 . 120021 ) . it is now 
generally accepted that the conditions for fragmentation 
are difficult to achieve i n the inner planet forming re- 
gions of protostellar disc s ([M atzncr & Levin 2005; R afikovl 
20051: iBolev et all 120061: IWhitworth fc Stamatellosl 120061 : 
Stamatellos fc Whitwortbll2008l ; iForgan et al .1120091) . A self- 
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gravitating protostellar disc is more likely to settle into a 
quasi-steady state in which the i nstability acts t o transport 
angu l ar momentum outwards (iGammid 2001; Rice et al 
l2003l ; lLodato fc Rice! [20041; IVorobvov fc Basul 120071) . It has 
been shown l|Lodato fc Riceil2004 120051 ) that in this quasi- 
steady state, the transport can be approximated as a local 
viscous process for all but the most massive discs. Since 
the disc will also be in thermal equilibrium, the strength 
of the angular mom entum transpo r t is also set by t he cool- 
ing rate of the dis c (lGammidl200ll; IRice et al-lliooll ). Using 
the above results, IRice fc Armitaga ( 20091 ) have shown that 
if self-gravity is the dominant transport mechanism, pro- 
tostellar discs will settle into quasi-steady states that are 
largely independent of the initial conditions. The resulting 
surface density profiles are quite steep wi th ~ 80 % of the 
mass w ithin 50 au located inside 10—20 au. lRice fc Armitage! 
(2009) also found that the quasi-steady mass accretion rates 
depend strongly on the disc mass and on the mass of the cen- 
tral star. Their simulations also suggested that no regions of 
the disc (which extended to 50 au) were susceptible to frag- 
mentation with the inner 10 — 20 au being pa rticularly sta- 
ble, consisten t with other recent calculations (Clarke 2009; 
lRafikov|[2009l ). 

The simulations carried out bv lRice fc Armitage! ([2009) 
assumed initial power-law surface density profiles that were 
evolved until a quasi-steady state - which we define as a 
state in which the mass accretion rate is approximately in- 
dependent of radius - was achieved. Although this gives 
information about the quasi-steady nature of such sys- 
tems, it doesn't tell us if such a state can be achieved 
under realistic conditio ns. Here, we extend the work of 
IRice fc Arniita gc (2009) by considering the evolution of self- 
gravitating, circumstellar discs formed through infall from 
a molecular cloud core. We find that these disc do settle 
rapidly into qu asi-steady states, with p roperties similar to 
that found by IRice fc Armitagd |2009). The mass accre- 
tion rates are also ge nerally consistent with observations 
l|Muzerolle et al.ll2005l ). although at early times the simu- 
lated accretion rates are somewhat higher than observed. 
Although the disc masses are generally higher than ob- 
served, a significant fraction of the mass is located in the 
optically thick inner regions and would not be easily de- 
tected using current tech niques. These results are consis- 
tent with the suggestion dArmitage. Livio fc Pri nglc 2001; 
IZhu. Hartmann fc Gammid 120091 ) that mass will pile-up in 
the inner regions of the disc and will drain episodically 
onto the star - potentially e xplaining FU Orionis outbursts 
l|Hartmann fc Kenvonl Il996l ) - when the temperature be- 
comes sufficiently high for MRI to operate. 

The paper is organised as follows. In Section 2, we de- 
scribe the basic model that builds the disc from an infalling 
molecular cloud core and evolves the disc self-consistently 
through, primarily, self-gravity. In Section 3, we describe 
the results, and in Section 4 we summarize our conclusions. 



2 BASIC MODEL 
2.1 Infall 



We b ase the model here on that of iLin fc Pringld 
(1990) which itself is based on the earlier work of 



ICassen fc Moosmanl (|l98ll ). We consider a molecular cloud 
core of mass M c , radius R c and angular velocity f2 c that is 
assumed to have an initially uniform density, p c . The core is 
allowed to collapse under gravity, with the collapse occuring 
on a characterist i c free -fall time of tft = (37r/32Gp c ) 1 ^ 2 . As 
m ILin fc Pringlel (ll990D. we assume that the collapse is ini- 
tially steady (|Larsonll 19691 ) and then, after one free-fall time, 
decreases smoothly to zero over a furth er free-fall time. The 
form of the accretion rate is therefore (ILin fc Pringld [19901 ) 



it/,. 



M(t) = { Iff [l + cos ^M] 




t < i ff 
iff < t ^ 2ta 
2t a sC t 



(1) 



which implies that the total amount of mass accreted is ac- 
tually 1.5M C . The infalling material is assumed to conserve 
angular momentum and consequently produces a centrifu- 
gally supported disc which, in the absence of angular mo- 
mentum transport, would have a maximum radius of a ma x = 
Q. 2 c Rt/GM c . We model the infall l|Cassen fc Moosman|[l98ll ) 
as the collapse of shells with each part of a shell striking the 
disc plane at the same time. The fraction of he mass of a 
shell, with initial radius r and enc losed mass M(< r ) , that 
strikes the disc within radius r is (|Lin fc Pringlel Yl 990) 



g(r) = 1 



where 



1-1- 

do 



1/2 



< r < a a , 



(2) 



(3) 



GM(< r ) ' 

At t = the initial timestep, dt, determines, through 
equation (fTJ, the amount of mass added and hence the 
amount of mass in the current shell. Since the core is as- 
sumed to have a uniform density, p c , this determines the 
radial extent of the initial shell. Equation @, together with 
equation (|3}, then determines how the mass in the shell is 
distributed in the disc. In subsequent timesteps, the radial 
extent of the next shell is then determined by ensuring that 
the amount of mass in the shell matches the amount of mass 
being added in that timestep, and the mass in that shell is 
distributed in the same way as in the first timestep. 

2.2 Disc evolution 

In order for the disc to evolve there needs to be some mech- 
anism for transporting angular momentum outwards, allow- 
ing mass to accrete onto the central star. Generally this is 
assumed to take the form of some kind of kinematic viscos- 
ity, v, which would cause an axisymmetric disc of surface 
density i) to evolve as (|Lvnden-Bell fc Pringlel Il974l ; 

|Pringldll98ll ) 



as 
dt 



3d_ 

r dr 



,1/2 



d_ 

dr 



(4) 



In the work of lLin fc Pringld j|l990h the kinematic vis- 
cosity is assumed to be the sum of two components. The 
first component is a standard "turbulent" viscosity given by 



(•») 



where c s is the disc sound speed, Q is the angular frequency, 
and a ss << 1 is a parameter that determines the efficiency 
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of th e angular momentum transport (Sha kura fe Sunvaevl 
Il973h . In astrophysical discs, the mechanism that provides 
this kinematic viscosity is generally accepted to be magneto- 
hydrodynamic (MHD) turbulence initiated by the magneto- 
rotational instability (jBalbus fe Hawlevlll99ll ). The second 
component, v s , was only non-zero if the disc was self- 
gravitating. 

An accret ion disc can become gravitational unstable if 
<|Toomrelll9"o4h 



Q 



ttGE 



(6) 



where k is the epicyclic frequency and can be replaced by 
the angular frequen c y, O, in a Keplerian disc. Following 
iLin fe Pringiel (ll987l ). lLin fe Pringld (|l990T ) assumed that u s 
was zero unless Q ^ Q c , where Q c was a critical value below 
which the disc became gravitationally unstable, and that v s 
increased as Q decreased. The evolution of self-gravitating 
discs h as, however, been studied in quite some detail re- 
cently lGammiell200ll;lRice et al.ll200d ; lLodato fe Ricdl2004l; 
iMeiia et al.ll2005l ; iDurisen et al. l2007t ) and it is now gener- 
ally thought that such discs will tend to settle into a state 
with constant Q and in which the heating and cooling rates 

balance. 

This work therefore differs from ILin fe Pringiel (Il990l ) 
in two ways. We firstly assume - at least during the ear- 
liest phases of star formation - that in general the disc is 
too cold and dense to have even the smal l degree of ion- 
izatio n needed to sustain MHD turbulence (|Blaes fc Balbud 
1 19941 ). In this case a BB = and disc self-gravity is the 
dominant mechanism for transporting angular momentum. 
We secondly assume that the disc will quickly settle into 
quasi-steady state with Q ~ 2. In this state the disc is as- 
sumed to be in thermal equilibrium with the radiative cool- 
ing balanced by heating occuring primarily through com- 
pressions and shocks. If we assume that the energy is dis- 
sipated locally, which appe a rs to be the case wh en Q ~ 1 
jBalbus fc Papaloizou|[l999l ; lLodato fc Rice] |2004| ) . the dis- 
sipation is equivalent to that occurring through an effective 
gravitational viscosity. Matching the dissipation rate with 
the known cooling rate allows the effectiv e gravitational a 
to be determined l|Rice fc Armitaed 120091 ) and the disc is 
then evolved using equation (|4]). 



the condition that Q — 2. We have the additional constraint 
that T c cannot be less than 10 K. We approximate the op- 
tical depth using 



dzn(p z ,T z )p z » Hk(p,T)p, 



(9) 



where k is the Rosseland mean opacity and p and T are 
an average density and temperature for which we use p — 
E/(2i?) and T = T c . The Rosseland mean opaci ties are de- 
termin ed using the analytic approximations from lBell fc Linl 
(1994). Once the optical depth is k nown, the cool ing func- 
tion, A, can be approximated using l|Hubenvlll990l ) 



3 K c oj l + r 2 ' 



(10) 



where T = 10 K is an assu med minimum temperat ure set by 
some background sources (|Stamatellos et al"1l2007l ). and the 
last term is introduced to smoothly int erpolate between op- 
ticall y thick and optically thin regions (jjohnson fc Gammid 
2003). The effective gravitational viscosity can then be de- 
termined by equating equations (0 and (j 10p . If, however, 
we rewrite the gravitational viscosity as v g — a g c? a /Q, and 
use that t coo \ — U /A, where U is the internal energy per 
unit area given by 



U 



(11) 



(12) 



7(7-1)' 
we get (|Gammidl200ll ) 

4 

Qfl ~ 9 7 (7-l)t cool fi" 

Although we are assuming that the primary transport 
mechanism is self-gravity, if regions of the disc were hot 
enough to be partially ionised, MRI may play a role in these 
regions. We therefore assume that if the temperature ex- 
ceeds ~ 1400 K MRI will operate and set a = a ss + a K = 



0.01. It has been suggested ( 


Armitage. Livio & Pringldl200ll; 


IZhu. Hartmann & Gammie 


2009) that this may explain FU 



Orionis outbursts. Self-gravity will transport mass to the 
inner disc where it will pile-up and cause the temperature 
in the inner disc to rise in order to remain gravitationally 
stable. Once hot enough for there to be partial ionisation, 
MRI will operate, rapidly draining the inner disc producing 
an FU Orionis-like outburst event. 



2.3 Determining the effective gravitational 
viscosity 

If we assume that a quasi-steady, self-gravitating disc intro- 
duces an effective gravitational viscosity, v g , the diss ipation 
in th e disc, D(R), per unit area per unit time is l|Bell fc Linl 
1 19941 ) 

D(R) = \v CJ Y.Q. 2 . (7) 

To determine the cooling rate we need to know the midplane 
temperature and the optical depth. The midplane tempera- 
ture can be determined using 

T c = (8) 
7«b 

where p = 2.4 is the mean molecular weight, m p is the 
proton mass, 7 = 1.4 is the specific heat ratio, fee is Boltz- 
mann's constant, and c s is the sound speed determined from 



2.4 Numerical method 

We consider an initial core mass, M c , with angular veloc- 
ity, n c and choose initially a small timestep, dt. The disc 
is modelled by a grid with 500 logarithmically spaced grid 
points extending from n n = 0.25 au to r out = 10a o . This 
does mean that the grid spacing is not the same for dif- 
ferent values of fi c , but does ensure that the outer edge of 
the disc should never reach the outer edge of the grid. The 
initial infall is modelled as described in Section 2.1 and we 
assume that any material that falls inside the inner edge of 
the disc falls directly onto the protostar. After the first small 
timestep we then have a small amount of mass on the proto- 
star, and the rest of the material distributed throughout the 
disc as described by equation ([2]). This means we now know 
the disc surface density, E. The midplane temperature, T c , 
is then determined using Q = 2 with the additional con- 
straint that T c > 10 K. The optical depth, r, cooling time, 
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t coo i, and effective gravitational visocsity, v g — a g c 2 /fl, can 
then be determined and the disc is evolved using equation 
. The amount of mass accreting onto the star is deter- 
mined from the difference between the disc mass before 
and after viscous evolution. The new mass distribution in 
the disc also changes the gravitational potential resulting 
in a lack of force balance in the disc. To achieve force bal- 
an ce material is m o ved in wards using the method described 
in iBath fe Pringld ljl98ll ). The infall is then repeated with 
again some mass falling directly onto the protostar and the 
rest being distributed throughout the disc. However, when 
adding mass to the disc, the added material will generally 
have a lower specific angular momentum than the mate- 
rial already in the disc. This gives rise to an inflow which 
we take into account u sing, again, the method described by 
IBath fc Priiiglel l|l98ll ). This once again changes the gravi- 
tational potential and so the lack of force balance must also 
be accounted for as before. The new disc surface density can 
then be used to recalculate the new gravitational viscosity 
and the process is repeated. 

After the initial timestep (in which we essentially pop- 
ulate the disc with a small amount of mass) , the subsequent 
timesteps are determined after the viscosity has been calcu- 
lated and are set to be a small fraction of the smallest diffu- 
sion time across a grid cell. We use a zero-torque boundary, 
S(ri n ) = 0, at the inner edge of the grid and although the 
outer boundary has no influence on the results, we set v r = 
at this location. 



3 RESULTS 

We consider core masses from 0.25 Mq to 5 Mq. In all cases 
we choose the radius of the cloud such that the resulting 
uniform density of p — 6 x 10~ 19 g cm -3 gives a free-fall 
time, tff, of 86900 years. For M c = 1 Mq this corresponds 
to a radius of 0.03 pc. The angular velocity is parametrised 
using / = fl c / \/Gp c and we vary / from / = 0.05 to / = 1.3, 
corresponding to angular velocities between 1.6 x 10 -15 rad 
s _1 and 2.6 x 10~ 13 rad s _1 . In many situations, the angular 
velocity is represented as the ratio of the total kinetic energy 
of rotation to the absolute value of the gravitational poten- 
tial energy, often referred to as /3. The relationship between 
f3 and / is /3 = f 2 /{2n) and so our chosen rotation rates 
correspond to j3 values between f3 = 0.0004 and /3 = 0.27. 



3.1 Quasi-steady evolution 

It has been suggested l|Rice fc Armitage! lioogh that a self- 
gravitating protoplanetary disc will settle into a quasi- 
steady state with a mass transfer rate that is approximately 
the same at all radii. iRice fc Armi tage (2009]), however, as- 
sumed initial power-law surface density profiles which were 
evolved until a quasi-steady state was achieved. Here, we 
self-consistently build and evolve the disc through the col- 
lapse of a uniform density molecular cloud core. Figure \T\ 
shows the disc properties at t = tg, t = 2tg, and t = 5fg 
for M c = 1 Mq and for / = 0.1. Each panel shows the disc 
surface density (solid line), central temperature (dash-dot- 
dot-dot line), effective gravitational a (dash-dot line), and 
mass accretion rate (dashed line) all plotted against radius 



from 0.25 au to 200 au. The left-hand y-axis is scaled for sur- 
face density and temperature, while the right-hand y-axis is 
scaled for a g and mass accretion rat e. The mass a ccretion 
rate is determined using M = 37r;/S l|Pringlelll98lh . Figure 
[2]is the same as Figure [1] except that / = 0.8 (i.e., the initial 
core rotation rate is faster in Figure [5] than in Figure (TJ. 

As suggested by I Rice fc Armitage! ((2009), the disc 
quickly settles into a quasi-steady state with a mass trans- 
fer rate that, apart from the very inner disc, is similar at 
all radii. The low mass transfer rate in the inner disc is due 
to the assumption that MRI ope rates when T c > 1400 K 
l|Zhu. Hartmann fc Gammiel [2009). which we approximate 
by setting a = a ss + a s = 0.01 when this condition is satis- 
fied. The disc therefore has short, episodic periods of rapid 
accretion (not shown in Figure [T] or Figure [2| draining mass 
from the inner disc, and much longer quiescent periods while 
the inner disc is replenished. We don't, however, claim that 
we are modelling the processes in the inner disc particularly 
accura tely. A mu c h mor e detailed analysis has been carried 
out bv lZhu et all (|2009h . 

For / = 0.1 (Figure [T] the disc is quite compact due to 
the slow rotation producing a small a max . The disc does ex- 
pand slightly with time, reaching an outer radius of ~ 100 
au at t = 5tf[. For / = 0.8 (Figure [2J the disc radius is 
much larger and actually extends initially beyond the 200 
au shown in Figure [2] In both cases, the d isc properties are, 
howev er, essentially the same as shown bv lRice fc Armitage! 
(2009) and remain largely unchanged over the 5 free-fall 
times shown in both Figure [1] and Figure [2] In this quasi- 
steady state, the disc has a reasonably steep surface density 
profile (E oc r -1,5 ). The a values, on the other hand can 
increase dramatically with radius due to the optically thick 
inner disc resulting in a very long cooling times - and hence 
small q values - and the optically thin outer disc resulting 
in short cooling times and large a values. 

The surface density profile is also not a pure power-law, 
and has a significant fraction of the mass in the inner disc 
(r < 10 au). The form of the surface density profile can 
be largely understood from the temperature profile. In the 
outer disc, the temperature is limited by our assumed mini- 
mum of 10 K. The temperature rises with decreasing radius, 
eventually reaching the temperature (T c = 170 K) above 
which ice mantles can no longer exist on the solid grains, 
generally known as the ice- or snowline. This produces an 
abrupt change in the opacity, the local cooling time, and 
consequently the effective gravitational viscosity. The ten- 
dency to evolve to a state with roughly constant M then 
produces a corresponding change in the surface density pro- 
file. The temperature then continues to rise with decreasing 
radius, eventually reaching 1400 K and becoming sufficiently 
ionised for MRI to operate. This sets an effective maximum 
temperature and also produces a subtle change in the sur- 
face density profile. A detailed description of how the surface 
density and temperarture profiles vary in the different opac- 
ity regimes is included in IClarka |2009). 

Not only does the quasi-steady nature of the disc persist 
for a number of free-fall times, it is also qualitatively the 
same for all our chosen parameters. Figure [3] shows the disc 
properties (as in Figure[T]and Figure[2]) at t — tff for / = 0.8 
and for M c = 0.5 M (top panel) and M c = 5 M (bottom 
panel). The form of the profiles is essentially the same as 
in Figures [1] and [2] A few quantitative differences are that 
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Figure 1. Disc properties for M cor c = 1 Mq for / = 0.1 at t = fff , 
t = 2tg, t = 5tff . Each panel shows the surface density S, (solid 
line), temperature (dash-dot-dot-dot line), effective gravitational 
a (dash-dot line), and mass accretion rate (dashed line). The left- 
hand y-axis is for S and temperature, while the right-hand y-axis 
is for a and mass accretion rate. This figure illustrates that the 
disc quickly settles into a quasi-steady state that persists for many 
free-fall times. 



Figure 2. The same as in Figure[T]except for / = 0.8. This again 
shows that the disc settles into a quasi-steady state that persists 
for many free-fall times. 



the surface density, the radius of the ice/snowline, and the 
radial extent of the region over which MRI can operate all 
increase with increasing core mass. Figure [2] also illustrates 
how the snowline moves in with time starting at ~ 50 au 
when t = tff and moving in towards 10 au at t = 5tg. 
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Figure 3. Disc properties for / = 0.8 and for M coro = 0.5 Mq 
(top panel) and M coro = 5 Mq. This figure illustrates that the 
quasi-steady nature of self-gravitating discs is largely the same for 
all core masses. There are some quantitative differences such as 
the surface density, the radius of the ice/snowline, and the radial 
extent over which MRI can operate all increase with increasing 
core mass. 



3.2 Mass accretion rates 

We compute the accretion rate onto the protostar by simply 
determining how its mass changes with time. The disc accre- 
tion rate is computed using M = SirvY, measured at a fixed 
location in the disc (generally 3 au for slowly rotating cores 
and 10 au for rapidly rotating cores) . Figure [4] shows how 
the disc accretion rate (dashed-line) and protostar accretion 
rate vary with time for M c — 1 Mq and for / = 0.1 (top 
panel) and / = 0.8 (bottom panel). In both cases shown the 
protostar accretion rate rises rapidly with time and peaks at 
t ~ tg (~ 10 5 years). It then drops significantly over the next 
freefall time and then changes more slowly from t = 2tg to 
t — 1 Myr. The protostar accretion rate also becomes quite 
variable after t = 2tg due to the accretion onto the proto- 
star occuring episodically when T c > 1400 K. In fact the 
protostar accretion rate has been averged over O.ltg and so 
the variability that our model produces is somewhat more 
extreme than shown in Figure [4] although we don't claim to 
be modelling this variability particularly accurately. 



The disc accretion rate (dashed line) is much smoother 
than the protostar accretion rate and, in Figure [4] has not 
been averaged. This illustrates that the disc is able to reach 
a quasi-steady state on timescales much shorter than the 
freefall time. What Figure [4] also shows is that for rapidly ro- 
tating cores (large /) the accretion onto the core is governed 
by accretion through the disc, while for slowly rotating cores 
the protostar grows initially through direct infall and only 
after infall ceases (t > 2tg) does disc accretion determine 
the rate at which mass is accreted onto the protostar. The 
accretion rate at t — 1 Myr also depends on the core rota- 
tion and increases with increasing /. From Figures [TJ [5] and 
[3] it is clear that the surface density profiles do not vary sig- 
nificantly with / or with time. What is different is the mass 
of the central protostar. The protostar mass increases with 
decreasing / and consequently (since we assume Q — 2) pro- 
duces disc temperatures that also decrease with decreasing 
/. This means that the cooling rate, effective gravitational 
viscosity, and consequently the disc accretion rate all de- 
crease with decreasing /. In fact for very slow rotating cores 
(/ < 0.08) it is difficult to sustain disc accretion beyond a 
few freefall times. 

Since we have carried out a large number of simu- 
lations with various core masses and with various rota- 
tion rates, we can consider how the mass accretion rate 
varies with protostar mass during the first 1 Myr. Figure [5] 
shows the mass accretion rate against protostar mass from 
our simulations (solid dots) compared wi th observed ac- 
cretio n rates for T Tau r i stars t aken from Gullbring et al.l 
(1 19981). I White fc Ghezl l|200l|) , ICalvet et all (|2004h . and 
iNatta. Testi fc Randichl (200€). The data points from our 
simulations are at t = tg (the peak accretion rate) and then 
at 50 randomly chosen times between t = tg and t = 1 Myr. 
The reasonably sparsely populated upper region illustrates 
how the accretion rate drops, within about a freefall time, to 
a value an order of magnitude or more lower than the peak 
value. The lower region is more densely populated, making 
it more likely that observed accretion rates would fall in this 
band. This likelihood is increased further due to the system 
being more heavily embedded during the first few freefall 
times, making it harder to measure accurate accretion rates 
at these early times. 

The accretion rates from the simulations don't compare 
particularly well with the observed accretion rates (open 
circles), with the observed accretion rates tending to be 
lower than the those from the simulations. The observed 
systems are, however, probably at a later stage of evolution 
than the simulated systems. The highest observed accretion 
rates are, however, consistent with what we would expect 
for the youngest systems, suggesting that self-gravitating 
transport can dominate until the earliest stages of the 
TTauri phase. It has also been suggested that the accre- 
tion rate should vary approximately with mass as M oc M 2 
IIMuzerolle et al. I l2005l ; INatta. Testi fc Randichl 120061 1 . The 
upper boundary of our simulated accretion rates, however, 
has a mass dependence closer to M oc M (illustrated by the 
sol id line in Figured] ) and i s consistent with that suggested 
bv lVorobvov fc Basuf (|2009^ and with the upper envelope of 
observed accretion rates (Hartmann et al . 2006). That the 
observed accretion rates are for protostars of various ages 
will introduce a scatter that could result in what appears to 
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Figure 4. Figure showing the accretion rate onto the protostar 
(solid line) and through the disc (dashed line) for M core = 1 
Mq and / = 0.1 (top panel) and / = 0.8 (bottom panel). The 
accretion rate initially increases with time, peaking at t ~ tg . The 
accretion onto the protostar is also quite variable and in fact has 
been averaged over O.ltg so is actually more variable than shown. 
The disc accretion is, however, much steadier and illustrates the 
quasi-steady nature of the disc evolution. The figure also shows 
that for high core rotation rates (/ = 0.8) accretion onto the 
protostar is governed by accretion through the disc, while for 
slow core rotation rates (/ = 0.1) the protostar grows primarily 
through direct infall during the first 2 free-fall times. 

be a steeper dependence on protostar mass than is actually 
the case. 

3.3 Disc fragmentation 

It has been proposed (Boss 1998, 2002) that discs that are 
sufficiently gravitationally unstable could fragment to pro- 
duce bound objects and that these bound objects could 
subsequently contract to form either gas giant planets 
or brown dwarf stars. What is now well known is that 
disc fragmentati on requires rapid cooling ( G ammiel l200ll ; 
iRice et all 120031 ) and that there is a minimum cooling 
time - which may depend on the actual cooling func- 
tion (Cos sins. Lodato fc Clarke] [20 09) - for which a self- 
gravitating disc can remain in a quasi-steady state with- 



Figure 5. Protostar accretions rates plotted against protostar 
mass from all our simulations (solid dots). The data points are 
at t = tff and at 50 random chosen times between t = tg and 
1 Myr. The sparsely populated region at the top illustrates how 
the accretion rate drops rapidly between t = iff and i = 2tff and 
suggests that we are not particularly likely to observe systems 
in this state. The open c ircles show observed ac c retion rates for 
TTau ri stars taken from iGullbring et all dl998h. IWhite fc Glie3 
feOOlh . lCaTvet et alj d2004h . iNatta. Testi fc Randichl (120061) . Al- 
though these don't compare particularly well with our calculated 
accretion rates, the highest observed accretion rates are consistent 
with what we would expect for the youngest systems, and many 
of the observed systems are probably at a later stage of evolution 
than the simulated systems. The solid line shows that the peak of 
the simulated accretion rates depend linearly on protostar mass 
(M oc M) and has the sa me relation as the uppe r envelope of the 
observed accretion rates l lHartmann et al. 2006). 

out fragmenting. The relationship between effective gravi- 
tational a s and cooling time (see equation (|12[0 allows us 
to define the fragmentation b oundary in terms of a g . It has 
been shown l|Rice et al.ll2005l ) that fragmentation occurs for 
a g > 0.06 and that this boundary is independent of the 
specific heat ratio 7. 

Figures [1] [2] and [3] show that in all cases, the a g values 
in the inner disc are well below that required for fragmen- 
tation (a g < 10~ 2 ), but that the outer parts of some discs 
could be susceptible to fragmentation. The requirement for 
fragmentation is that Q ~ 1 and a g > 0.06. In our sim- 
ulations we assume Q settles to Q = 2, so we don't quite 
satisfy the first condition, but we can at least identify the 
regions where the a s condition is satisified and therefore 
where Q will decrease and fragmentation is most likely. Fig- 
ure [6] shows the range of radii where fragmentation could 
occur at t = tfi (triangles) and t = 1.5ts (diamonds), for 
various / values. Each vertical line is for a single core mass 
and we consider core masses of M c = 0.25 Mq (solid line), 
M c = 0.5 M Q (dashed line), M c = 1 M (dash-dot line), 
M c = 2 M Q (dash-dot-dot-dot line), and M c = 5 M Q (long 
dash line). 

Figure [6] shows that only the outer disc (beyond ~ 40 
au) has conditions suitable for fragmentation, cons istent 
with numerical simulations ( Sta matellos et all 120071) an d 
possibly also with some observations jGreaves et alj|2008h . 
The inner radius of the fragmentation zone increases with 
core mass and for a core mass of 1 Mq, which produces a 
central protostar with a mass close to 1 Mq , is 60 — 70 au, 
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occurs if £/3 < 0.12. If we assume an initially i sother mal 
core, with sound speed c s , then (see lWalch et al.l (|2009l )) 



Figure 6. Figure showing the range of radii where fragmentation 
could occur (defined as the region where a > 0.1). We find that 
the fragmentation conditions are never satisfied when / ^ 0.1 and 
are always satisfied when / 0.25. To illustrate that the fragmen- 
tation conditions can be satisfied for many dynamical times, the 
fragmentation range is plotted at t = tg (triangles) and t = 1.5% 
(diamonds). The core masses considered are Mcore = 0.25 Mq 
(solid line), M corc = 0.5 M Q (dashed line), M core = 1 M (dash- 
dot line), M core = 2 Mq (dash-dot-dot-dot line), and M corc = 5 
Mq (long-dash line). The inner fragmentation radius increases 
with core mass from ~ 40 au for Mcore = 0.25 Mq to ~ 100 
au for M core = 5 Mq. The radial range also increases with in- 
creasing core mass, but this should be interpreted cautiously as 
the outer disc could be susceptible to ionisation fro m cosmic rays 
that may stimulate MRI turbulence jRafikovll20"09h and stabilise 
these regions against fragmentation. 



consistent with IClarkd (|2009l ). The radial range that can 
undergo fragmentation is quite small for low core masses 
(10 — 20 au for M c — 0.25 Mq) but increases dramatically 
with core mass. This result needs to be interpreted some- 
what cautiously as the outer parts of these discs could be 
susceptible to ionizatio n through cosm ic rays that may stim- 
ulate MRI turbulence (|Rafikovi r2009) and stabilise these re- 
gions against fragmentation. 

Although we can't really predict the possible masses 
of the objects that may f orm via fragmentation, 
iKratter. Murrav-Clav fc Youdirj (|2009) suggest that the ini- 
tial fragmentation masses will be ~ 1 Mj up itcr- The large 
amount of mass in these outer regions, however, suggests 
that such fragments will continue to grow and that the iso- 
lation mass can easily be > 100 Mj up iter- Fragmentation in 
the outer parts of these discs is therefore unlikely to form 
planetary mass objects and will probably form sub-stellar 
or stellar mass companions. 

None of our simulations with / = 0.1 or less had 
any regions susceptible to fragmentation. This was essen- 
tially due to these slow rotation rates producing compact 
discs with most of the mass in the optically thick inner 
disc and very little mas s in the optically thin outer disc . 
It has been suggested (iMivama. Havashi fc Narital I1984T ) 
that fragmentation depends on both the ratio of the ther- 
mal to gravitational energy £ and the ratio of the en- 
ergy of rotation to the gravitatio nal energy [3 = f 2 /(2n). 
IMivama. Havashi fc Narital (| 1984T ) find that fragmentation 
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Since we assume all our cores have the same density and 
hence freefall time, the £ values of our cores vary with core 
mass. Since fragmentation is possible for all core masses with 
/ ^ 0.25 and for none of the core masses when / ^ 0.1, 
this suggests that what primarily determines if a disc will 
fragment is the initial rotation rate of the cloud core. The 
reason for this is that, in our simulations, the discs very 
quickly settle into a quasi-steady state, with Q — 2, that 
is in thermal equilibrium. This quasi-steady state is essen- 
tially independent of the initial thermal properties of the 
core and depends largely on how much mass is in the sys- 
tem. If the initial rotation is such that some of this mass can 
be deposited in the outer, optically thin parts of the disc, 
then fragmentation is quite likely. Figure [6] also shows that 
the fragmentation region remains largely unchanged for at 
least 0.5tff (from t = ts to t — 1.5tg). We find, in fact, that 
conditions suitable for fragmentation actually persist for a 
freefall time or longer (which is many dynamical times even 
in the outer disc), but are generally no longer satisfied after 
~ 2 freefall times for the lower mass cores and after 3 — 4 
freefall times for the higher mass cores. 



3.4 Disc and protostar masses 

Figure[7]shows the evolution of the disc and protostar masses 
against time for M core = 1 Mq and for / = 0.1 (top panel), 
/ = 0.25 (middle panel), and / = 0.8 (bottom panel). In all 
cases, the solid line shows the evolution of the core mass, 
while the dashed lines show the disc mass within 10 au, 100 
au, and the total disc mass. For / = 0.1 the disc does not 
extend beyond 100 au and so the mass inside 100 au and 
the total disc mass are the same. The figure illustrates that 
as the core rotation rate increases, the core mass decreases 
and more and more mass is deposited at large radii (r > 100 
au). What these figures also illustrate is that the mass inside 
10 au and the mass inside 100 is only weakly dependant on 
the rotation rate, again illustrating the quasi-steady nature 
of these discs. 

Figure [8] is the same as Figure [7] except it is for / = 0.8 
and Mcorc = 0.5 M (top panel), M cola = 2 Mq (middle 
panel) and Af CO rc = 5 Mq (bottom panel). It is clear from 
both Figure [7| and Figure [8] that in all cases, there is an 
epoch when the disc mass exceeds the mass of the central 
object. For small / values, this epoch is quite short (< ltg) 
while for larger / values, the disc may remain more massive 
than the central star for many freefall times. The outer parts 
(r > 50— 100 au) are, however, susceptible to fragmentation 
and may break up to produce companions. The inner parts 
(r < 100 au) of these discs are only generally comparable in 
mass to the central protostar for t < 2tg. This does, however, 
suggest that our local approximation may not be strictly 
valid, and that the effective gravitational a will not be de- 
termined by the local cooling rate. I t has been sugg ested 
|Kratter. Matzner fc Krumhol j [ioosl ; IVorobvovll2lx)9r ) that 
there should be two self-gravitating a parametrisations, one 
when the disc-to-star mass ratio is small, and the other when 
the disc-to-star mass ratio is large. This may be the case, 
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but we assume that even if the disc is massive relative to 
the central star, the system would still settle into a state 
with constant Q and be in thermal equilibrium. We don't 
know if, in such a state, the local dissipation rate will be 
set by the local cooling rate, but we assume that such an 
assumption is at least reasonable at this stage. The figures 
also illustrate that in a quasi-steady, self-gravitating state at 
least 50 % of the disc mass within 100 au is located inside 10 
au. This large amount of mass in the inner disc could have 
implications for, and could aid, planet formation. Detecting 
such massive inner discs is, however, extremely difficult as 
they are very optically thick. 

As discussed in the previous section, it is also likely 
that for / ^ 0.25 the outer disc will be very unstable and 
may fragment to produce companions. The exact radius at 
which this occurs depends on the mass of the central star, 
but is typically between 60 and 100 au. A significant fraction 
of the material in the outer disc (r > 100 au) is therefore 
likely to be converted into companions which will truncate 
the outer disc at a radius of 50 — 75 au, depending on where 
the fragmentation occurs. A standard way to determine disc 
masses observationally is to fit the mass by modelling the 
spectral energy dist ribution (SED) (jRobitaille et al.l 120061 ; 
I Andrews et afll2009l ). The results from such modelling can, 
however, be ambiguous l|Eisner et al.l 12005) and, in partic- 
ular, a lot of mass can be hidden in optically thick regions 
of the disc. The enhanced mass in the optically thick inner 
10 au should therefore not contribute significantly to the 
SED, and the discs in Figures [7] and [8] would consequently 
appear to have a mass - determined through SED modelling 
- roughly equal to the difference between the mass within 
100 au and the mass inside 10 au. For TTauri-like protostar 
masses (e.g., Figure[7|, this would be 0. 1 — 0.2 Mr consiste nt 
with that obtained for class I sources (E isner et all [2005). 

If we only consider the disc mass within 100 au, since 
the outer disc is likely to be removed through fragmentation, 
Figure [8] and the bottom panel of Figure [7] also shows that 
the disc-to-protostar mass ratio increases with decreasing 
protostar mass, consistent with observations in the Orion 
Nebula cluster th at massive discs are more common around 
lower-mass stars (|Eisner et al.l 120081 ) . For core masses be- 
low ~ 1 Mq, however, the disc mass is very similar to the 
mass of the central star and again brings into question our 
assumption that the transport through self-gravity can be 
treated as a local process. The inner fragmentation radius 
for these lower-mass cores (see Figure [6} is, however, also 
closer to 50 au than 100 au potentially converting more of 
the mass in the disc into companions. The relatively low 
disc-to-star mass ratio for protostars with masses above a 
few solar masses may also help to explain why most discs 
around relatively massive stars disappear on much shorter 
timescales than discs around solar-like TTauri stars. 

3.5 Additional sources of viscosity 

The discussion in the section above suggests that fragmen- 
tation in the outer disc and the presence of an optically 
thick inner disc may result in very young protostars with 
TTauri-like masses appearing to have disc masses compara- 
ble to that observed. Figure [S] and Figure 2] however, show 
that at t = 1 Myr, the accretion rates around these pro- 
tostars is between 10 -8 and 10 -7 Mq yr _1 . At this rate, 
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Figure 7. Protostar (solid line) and disc (dashed lines) masses 
plotted against time for Mcore = 1 Mq and for / = 0.1 (top 
panel), / = 0.25 (middle panel), and / = 0.8 (bottom panel). 
The different dashed lines show the disc mass inside 10 au, 100 
au, and the total disc mass. The figure illustrates that as the 
core rotation rate increases, the core mass decreases and more 
and more mass is deposited at large radii in the disc. The figure 
shows that in all cases a significant fraction (~ 50 %) of the mass 
inside 100 au is located in the optically thick inner 10 au. 
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Figure 8. Protostar (solid line) and disc (dashed lines) masses 
plotted against time for / = 0.8 and for M CO rc = 0.5 Mq (top 



panel), M a 
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5 M Q . The 



different dashed lines show the disc mass inside 10 au, 100 au, and 
the total disc mass. The figure shows that a significant fraction 
(~ 50 %) of the mass inside 100 au is located in the optically 
thick inner 10 au. Considering only the disc mass inside 100 au 
(since the outer disc is potentially susceptible to fragmentation) 
the disc-to-protostar mass ratio decreases with increasing core 
mass which is consistent with observations in the Orion Neb- 
ula cluster ijEisner et al J 1200 8). and could help to explain why 
discs around relatively massive stars disappear on much shorter 
timescales than discs around lower-mass TTauri stars. 



removing the final disc material (Mdi sc > 0.3 Mq within 
100 au according to Figures [7] and [8]) would take in ex- 
cess of 10 7 yrs. This is significantly longer than the disc 
dispersal timescale of ~ 5 Myr suggested by observations 
of disc fractions in star fo r ming regions of d i fferent ages 
ijHaisch. Lada fc LadalboOll ). iRice fc Armitage! (|2009l ) have 
also shown that the mass accretion rate also drops dramat- 
ically with disc mass, becoming very small for disc masses 
less than a tenth of the protostar mass. If self-gravity is the 
only transport mechanism then the timescale for removing 
the disc material would be much longer than observed. This 
suggests that there should be an additional transport mech- 
anism. The most likely candidate is M RI operating in the 
upper layers of the disc (|Gammielll996t ). 



4 DISCUSSION AND CONCLUSIONS 

We consider here the evolution of protostellar discs that form 
through the collapse of molecular cloud cores and in which 
the primary transport mechanism is self-gravity. We assume 
that the discs settle quickly into a marginally gravitationally 
stable state with Q — 2 and are then in thermal equilibrium. 
In this state angular momentum transport is driven by an 
effective gravitational viscosity that is calculated using the 
assumption that energy is dissipated at a rate equal to the 
rate at which the disc los es energy through radiative cool - 
ing. Following the work o f Armitag e. Livio fc Pringlel (|200ll ) 
and lZhuTH artmann fc Gammid (|2009h we also assume that 
MRI l|Balbus fc Hawlevlll99lh will operate if the disc is hot 
enough (T > 1400 K) to be partially ionised. This occurs 
only in the inner disc and is episodic - draining material 
from the inner disc and then turning off until the inner disc 
has been replensished by material from the outer disc - and 
has therefore been pro posed as a mechanism for e xplaining 
FU Orionis outbursts (|Hartmann fc Kenvonlll996l ). 

Our simulations consider a range of cores masses (0.25 — 

5 Mq) and a range of rotation rates (0.05 ^ / ^ 1.3 with 
/ = fi c / y/Gpc)- The primary results are as follows 



• The discs settle quickly into a quasi-steady state with a 
reasonably steep surface density profile and with a lot (~ 50 
%) of the mass within 100 au located inside 10 — 20 au. 
The disc-to-star mass ratio decreases with increasing proto- 
star ma ss, consistent with observations in the Orion Nebula 
cluster (jEisner et al,l 120081 ) . Although the disc masses tend 
to be higher than observed, with so much mass located in 
the inner optically thick regions of the disc, most of this 
mass would be difficult to detect with current techniques. 

• Although the mass accretion rates are initially higher 
than those observed, these high accretion rates only per- 
sist for ~ 1 free-fall time and quickly drop to values sim- 
ilar to that observed. The simulations also suggest that 
the accretion rate varies less strongly with protostar mass 
(M oc Ma.) than sug g ested by observations ( M oc M» ) 
( Muzerolle et all 12003 ; iNatta. Testi fc Randicbl l200rj ) . but 
is consistent with the upper bound ary of the observed ac- 
cretion rates l|Hartmann et al]|2006l ). 

• In none of our simulations did the inner disc [r < 40 
au) have conditions suitable for fragmentation (a g > 0.06). 
In some cases, the outer disc was susceptible to fragmenta- 
tion, with the primary factor determining if a disc could 
fragment being the rotation rate of the molecular core. 
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The large amount of mass in these discs, however, sug- 
gests that fragmentation is more likely to result in sub- 
stellar or stellar companions that pla netary mass objects 
|Kratter. Murrav-Clav fc Youdir] [2009), None of our simu- 
lations with / ^ 0.1 satisfied the fragmentation conditions 
while it was satisifed for all those with / ^ 0.25. The in- 
ner radius of the fragmentation region increased from ~ 40 
au for Af corc = 0.25 M to ~ 100 au for M corc = 5 M 
as did the radial range over which fragmentation could oc- 
cur, although the outer regions of some of these discs could 
be stabilised a gainst fragmentation by cosmic ray ionisation 
l|Rankovll2009h . 

• The mass accretion rate depends strongly on the disc 
mass and in all cases drops below 10 -7 Mq yr _1 when the 
disc mass is still quite substantial (Afdi sc > 0.3 M©). If self- 
gravity were to remain the primary transport mechanism, 
disc clearing timescales would be significantly longer than 
that observed. This suggests that an additional tranport 
mech anism, such as MRI occuring in the upper layers of the 
disc (|Gammid[l996l ). must also operate. We propose that 
this additional mechanism is likely to be negligible initially 
and become more effective with time. This is probably the 
case for layered accretion which should become more effi- 
cient as the disc mass decreases and is consistent with the 
requirement that smalls dust grains must be depleted before 
MRI can operate effectively l|Sano et al. 2000; Sa no fc Stonel 
120021 ; lllgner fc Nelson|[2006l ). but the pile-up of material in 
the inner regions of the disc (ultimately resulting in episodic 
FU Orionis-like outbursts) would also not occur if the addi- 
tional transport mechanism were too efficient at early times. 

Ultimately it appears that transport driven by self- 
gravity can explain protostar formation and disc evolution 
at early times. Our simulations suggest that disc masses may 
- at early times - be higher than suggested by observations 
but these large disc masses also suggest that an additional 
transport mechanism must dominate at later times (t > 1 
Myr) to remove the remaining disc material within observed 
disc lifetimes. The pile-up of material in the inner regions 
of the disc - which may explain FU Orionis outburst - may 
also play an important role in the subsequent formation of 
planets. 



ACKNOWLEDGEMENTS 

P.J. A. acknowledges support from the NSF (AST- 
0807471), from NASA's Origins of Solar Systems pro- 
gram (NNX09AB90G), and from NASA's Astrophysics The- 
ory and Fundamental Physics program (NNX07AH08G). 
W.K.M.R. acknowledges support from the Scottish Universi- 
ties Physics Alliance (SUPA). The authors would also like to 
thank the Isaac Newton Institute for Mathematical Sciences 
for their hospitality during the Dynamics of Discs and Plan- 
ets Programme, and would like to acknowledge useful discus- 
sions with Lee Hartmann, Dick Durisen, Giuseppe Lodato 
and Neal Turner. 



REFERENCES 

Andrews S.M., Wilner D.J., Hughes A.M., Chunhua Q., 
Dullemomd CP., 2009, ApJ, 700, 1502 



Armitage P.J., Livio M., Pringle J.E., 2001, MNRAS, 324, 705 

Balbus S.A., Hawley J.F., 1991, ApJ, 376, 214 

Balbus S.A., Papaloizou J.C.B., 1999, ApJ, 521, 650 

Bath G., Pringle J.E., 1981, MNRAS, 194, 962 

Bell K.R., Lin D.N.C., 1994, ApJ, 427, 987 

Blaes O.M., Balbus S.A., 1994, ApJ, 421, 163 

Boley A.C., Meji'a A.C., Durisen R.H., Cai K., Pickett M.K., 

D'Alessio P., 2006, 651, 517 
Boss A.P., 1998, ApJ, 503, 923 
Boss A. P., 2002, ApJ, 576, 462 

Calvet N., Muzerolle J., Bricefio C, Fernandez J., Hartmann L., 
Saucedo J. L., Gordon K. D., 2004, AJ, 128, 1294 

Caselli P., Benson P.J., Myers P.C., Tafalla M., 2002, ApJ, 572, 
238 

Cassen P.M., Moosman A., 1981, Icarus, 48, 353 
Clarke C.J., 2009, MNRAS, 396, 1066 

Cossins P., Lodato G., Clarke C, 2009, MNRAS, in press 
Durisen R.H., Boss A. P., Mayer L., Nelson A.F., Quinn T., Rice 
W.K.M., 2007, in Reipurth B., Jewitt D., Keil K., eds, 
Protostars and Planets V. University of Arizona Press, 
Tucson, p701 

Eisner J. A., Hillenbrand L.A., Carpenter J.M., Wolf S., 2005, 
ApJ, 635, 396 

Eisner J. A., Plambeck R.L., Carpenter J.M., Corder S.A., Qi C, 

Wilner D., 2008, ApJ, 683, 304 
Forgan D., Rice K., Stamatellos D., Whitworth A., 2009, 

MNRAS, in press 
Gammie C.F., 1996, ApJ, 457, 355 
Gammie C.F., 2001, ApJ, 553, 174 

Greaves J.S., Richards A.M.S., Rice W.K.M., Muxlow T.W.B., 

2008, MNRAS, 391, L74 
Gullbring E., Hartmann L., Briceno C, Calvet N., 1998, ApJ, 

492, 323 

Haisch K.E., Lada E.A., Lada C.J., 2001, AJ, 121, 2065 
Hartmann L., D'Alessio P., Calvet N., Muzerolle J., 2006, ApJ, 
648, 484 

Hartmann L., Kenyon S.J., 1996, ARA&A, 34, 207 
Hubeny I., 1990, ApJ, 351, 632 
Ilgner M., Nelson R.P., 2006, A&A, 445, 205 
Johnson B.M., Gammie C.F., 2003, 597, 131 
Kratter K.M., Matzner CD., Krumholz M.R., 2008, ApJ, 681, 
375 

Kratter K.M., Murray-Clay R.A., Youdin A.N., 2009, ApJ, 

submitted 
Larson R.B., 1969, MNRAS, 145, 271 
Laughlin C, Bodenheimer P., 1994, ApJ, 436, 335 
Lin D.N.C, Pringle J.E., 1990, MNRAS, 358, 515 
Lin D.N.C, Pringle J.E., 1987, MNRAS, 225, 607 
Lodato C, Rice W.K.M., 2004, MNRAS, 351, 630 
Lodato C, Rice W.K.M., 2005, MNRAS, 358, 1489 
Lyndon-Bell D., Pringle J.E., 1974, MNRAs, 168, 603 
Matzner CD., Levin Y., 2005, ApJ, 628, 817 
Meji'a A.C, Durisen R.H., Pickett M.K., Cai K., 2006, ApJ, 619, 

1098 

Miyama S.M., Hayashi C, Narita S., 1984, ApJ, 279, 621 
Muzerolle J., Luhmann K.L., Briceno C, Hartmann L., Calvet 

N., 2005, ApJ, 625, 906 
Natta A., Testi L., Randich S., 2006, A&A, 452, 245 
Papaloizou J.C.B., Nelson R.P., 2003, MNRAS, 339, 983 
Pringle J.E., 1981, ARA&A, 19, 137 
Rafikov R.R., 2005, ApJ, 621, L69 
Rafikov R.R., 2009, ApJ, submitted. 
Rice W.K.M., Armitage P.J., 2009, MNRAS, 396, 2228 
Rice W.K.M., Armitage P.J., Bate M.R., Bonnell LA., 2003, 

MNRAS, 339, 1025 
Rice W.K.M., Lodato C, Armitage P.J., 2005, MNRAS, 364, 

L56 



12 W.K.M. Rice, J.H. Mayo & Philip J. Armitage 



Robitaille T.P., Whitney B.A., Indebetouw R., Wood, K., 

Denzmore P., 2006, ApJS, 167, 256 
Sano T., Miyama S., Umebayashi T., Nakaon T., 2000, ApJ, 

543, 486 

Sano T., Stone J., 2002, ApJ, 577, 534 

Shakura N.I., Sunyacv R.A., 1973, A&A, 24, 337 

Stamatellos D., Whitworth A.P., Bisbas T., Goodwin S., 2007, 

A&A, 475, 37 
Stamatellos D., Whitworth A.P., 2008, A&A, 480, 879 
Terebey S., Shu F.H., Cassen P., 1984, ApJ, 286, 529 
Terquem C.E. J.M.L. J., 2008, ApJ, 689, 532 
Toomrc A., 1964, ApJ, 139, 1217. 
van Boekel R., et al., 2004, Nature, 432, 479 
Vorobyov E.I., 2009, New Astron., 15, 24 
Vorobyov E.I., Basu S., 2009, ApJ, in press 
Vorobyov E.I., Basu S., 2007, MNRAS, 381, 1009 
Walch S., Burkert A., Whitworth A., Naab T., Gritschneder M., 

2009, MNRAS, in press. 
White R.J., Ghez A., 2001, ApJ, 556, 265 
Whitworth A.P., Stamatellos D., 2006, 458, 817 
Zhu Z., Hartmann L., Gammic C.F., 2009, ApJ, 694, 1045 
Zhu Z., Hartmann L., Gammie C, McKinney J.C., 2009, ApJ, 

701, 620 



